Renormalized mean-field theory for a two-component Fermi gas with s-wave 

interactions. 
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A method is introduced to renormalize the zero-range interaction for use in mean-field and many- 
body theory, starting from two-body calculations. The density-renormalized delta-function interac- 
tion is then applied using mean-field theory to a two-component fermion gas, and compared with 
diffusion Monte Carlo simulations and conventional mean-field calculations. In the unitarity limit, 
the equation of state exhibits the expected behavior fj, oc p 2//3 , with a parameter f3 = —0.492, which 
is consistent with recent experiments^, 000 • 

PACS numbers: 



I. INTRODUCTION 

The low energy scattering of fermionic atoms con- 
trols the structure and dynamics of an ultracold quan- 
tum Fermi gas. When the scattering length ao between 
fermions in different internal spin states is tunable, for in- 
stance, in an external magnetic field, it becomes possible 
to study the crossover between BCS-type superfluidity of 
momentum pairs and Bose-Einstein condensation (BEC) 
of molecules. 

In recent years, the BCS-BEC crossover problem has 
become experimentally accessible 0, IE 13, enabling 
sharp tests. The BCS theory has been successful in ex- 
plaining superfluidity in Fermi gases, but this theory is 
incomplete because it neglects the Hartree term of the 
interaction, Airap/m, where p is the density of one spin 
component. Comparatively little research has considered 
Fermi gases including the Hartree term, with the pri- 
mary regime studied being the perturbative case where 
pa^ << 1. This is usually referred as the "normal state" 
of the gas, in contrast with the superfluid state. Quan- 
tum Monte Carlo (QMC) simulations include both the 
Hartree term and pairing physics, but a complete theory 
that contains both ingredients is still required. 

When the range of interaction is much smaller than the 
inter-particle distance, the potential can be replaced by 
a delta function interaction although this must be done 
with caution because the delta function interaction is too 
singular to be exactly solvable, even in principle. In the 
weakly-interacting limit, pa^ « 1, the coupling param- 
eter in the delta function interaction is proportional to 
the two-body scattering length ao; this is known as the 
Fermi pseudo-potential 0. Using this approximation, 
mean field theories have been applied to Fermi gases 

nana. ^ use of this approximation in strong in- 
teracting regimes leads to an unphysical collapse of the 
Fermi gas. To go beyond the Fermi approximation for 
the pseudo-potential, it is crucial to renormalize the cou- 
pling constant; the purpose of this paper is to introduce a 
new and convenient way to achieve this renormalization. 

Full diagonalization of a Hamiltonian with delta- 
function interactions requires a momentum cut-off renor- 
malization even in the weak interacting limit. This type 



of renormalization has been carried o ut, for example, in 
Refs.[H EE O EE [IE E3, EE E3, ISO], just to name 
a few such studies. However, such a renormalization is 
unnecessary at small or modest scattering lengths, when 
treated by mean-field theories with zero-range interac- 
tions because they are well behaved in this limit. To go 
beyond the weakly-interacting limit of mean-field theory 
we propose a density-dependent renormalization of the 
coupling parameter that is intended to apply even in the 
long- wavelength limit. A density-dependent renormaliza- 
tion for a 2-component deg enerate Fermi gas has been re- 
cently proposed in Ref. [2l| to explain the stability of this 
system in the strong interacting limit; a functional form 
for the effective scattering length a e ff was designed to 
give the expected behavior (i.e. as determined by QMC 
and other calculations) in both the weak interaction limit 
and the unitarity limit. Our renormalization strategy is 
different, in that we present a method to calculate a e // 
by using the exact energies of two particles in a trap. 
We compare the many-particle predictions obtained us- 
ing the renormalized interaction potential with diffusion 
Monte Carlo and alternative mean-field calculations, and 
find that our renormalization automatically gives the cor- 
rect behavior in both the strong and weak interaction 
limits, for both positive and negative scattering lengths, 
without imposing this constraint at the outset. 

This paper is organized as follows. In Section II we 
develop the renormalization procedure and show that a 
simple 2-paramcter analytical formula can be utilized to 
an excellent approximation over the whole range from 
positive to negative two-body scattering lengths. Section 

III applies the renormalization to many-particle mean- 
field theory and presents some of its predictions. Section 

IV compares the results obtained using our renormalized 
mean-field theory with quantum Monte Carlo calcula- 
tions and with perturbative mean-field calculations. 



II. RENORMALIZATION PROCEDURE 

Through this paper we consider a system of equal mass 
fermions in a spherically symmetric harmonic oscillator 
trap at temperature T = 0. The Hamiltonian that we 
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adopt is 
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This Hamiltonian cannot be diagonalized exactly, since 
the delta function interaction is too singular and would 
produce divergent results. J22J The level of approximation 
we adopt to diagonalize the many-body Hamiltonian is 
the same we use to solve the corresponding Schrodinger 
equation for the two-body system. For example, if we 
want to diagonalize this Hamiltonian in the RPA approx- 
imation, we would use RPA for the two-body system and 
obtain the renormalization through the matching proce- 
dure explained later in this Section. An explicitly corre- 
lated wavefunction or an extensive configuration interac- 
tion (CI) wavefunction can produce divergent results and 
a momentum cutoff renormalization is necessary. Since 
we want to obtain and apply a density renormalization 
without the necessity of introducing a momentum renor- 
malization, we need to carefully select the level of approx- 
imation in the wavefunction. The Hartree-Fock (HF) 
approximation does not introduce explicit interparticle 
correlations, as the only correlations included are "ex- 
change correlations" from the Pauli exclusion principle. 
This makes the HF wavefunction a suitable approxima- 
tion to adopt in our renormalization technique, since it 
does not require a momentum renormalization. 

To obtain the renormalized scattering length we solve 
(|T|) for two opposite-spin fermions in the HF approxima- 
tion. The ground state energies of this approximation 
are matched with the exact energies of the system for 
different values of the bare two-body scattering length 
a . From this procedure we extract the functional de- 
pendence of a e /f on a,Q. The spectrum of two opposite- 
spin fermions in a trap having a specified scattering 
length an and zero-range interactions can be determined 
exactly p3. 123. l25l l26j] and (JTJ can be solved numerically 
for two particles using a Hartree-Fock wavefunction. 

Dimensional analysis suggests that, in an infinite, uni- 
form Fermi gas, where the range of the two-body interac- 
tion is much smaller than both the average interparticle 
distance and the bare scattering length an, the only pa- 
rameter that characterizes the behavior of the system is 
the dimensionless combination kja^ of the Fermi momen- 
tum and ao- Throughout this paper the Fermi momen- 
tum is defined as kt = {Sir 2 p) 1 / 3, where the density is 
just the one-spin component density. If we were apply- 
ing the renormalized scattering length a e f / to an infinite 
uniform system, the only relevant parameter would be 
kfa e ff. This suggests that kfa e ff has to be a function 
of fc/ao- So, we propose the following functional depen- 
dence, 



a eff 



CO/ao) 



(2) 



We will see below that the renormalization function, 
£(fc/ao), will have the desired behavior in the limiting 
cases, becoming independent of fc/ao in the unitarity 
limit (| do I — ► 00) and reproducing the relation a e ff = o-o 
in the weak interacting limit (fc/dn -C 1). We consider 
that Eq.J2J holds even with the inclusion of a trapping 
potential. The renormalized scattering length a e ff can 
be viewed as accounting, in some way, for the correlations 
neglected in the mean- field wavefunctions. 

In choosing to extract a e ff from a two-body system, 
we are implicitly assuming that two-body correlations are 
the most important in the many-body system. This as- 
sumption is reasonable for two-spin component fermions 
with short-range interactions because the probability of 
finding more than two fermions close enough to interact 
is usually negligible. 



A. Exact Energies 

This subsection reviews the exact results for two parti- 
cles in a trap interacting through a zero-range pseudopo- 
tential. Then, we will select from the energy spectrum, 
the energy branch that is of interest to study in the renor- 
malization procedure. 

Consider two particles of mass m interacting through a 
two-body potential V(r) in a spherically symmetric trap. 
If the effective range of the two-body potential is much 
smaller than the caracteristic length of the trap, the low- 
lying energy levels depend only on the scattering prop- 
erties of the potential and not on its shape. Under this 
condition, the two-body potential can be replaced by a 
pseudopotential of the form [2j| : 



v(r) = <H r )ir r ' 

m or 



(3) 



where a(E) is the energy-dependent scattering length. 
For ultra-cold gases, the energy dependence on the scat- 
tering length can be neglected, so a(E) can be replaced 
by its energy-independent limit an- 

Two particles in a trap with this pseu do p otential have 
been considered previously [23|, Hj, |2|| |26j . After sepa- 
rating the center of mass and relative coordinates, the 
problem reduces to two independent one-dimensional 
Schrodinger equations. The pseudopotential J3J is intro- 
duced as a boundary condition in the relative coordinate 
Schrodinger equation. The energies are 



E, 



CM 



E 



where Ecm = (ncM + 3/2)/kj and 
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a 



(4) 



(5) 



The trap length ato is defined in this case as d/j = 
y/h/mu. Equations Q and JSJ reproduce the com- 
plete spectrum of the system with zero relative angu- 
lar momentum. Figure ^ shows the spectrum of E re i 
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(Eq. (0) as a function of the scattering length. The 
lowest curve of the spectrum shown describes the forma- 
tion of a molecule, where interparticle correlations are 
fundamental. One anticipates that a HF wavefunction 
would be a terrible approximation to such a state in 
which the two atoms are bound together to form a bound 
molecular eigenstate. Since in this work we do not con- 
sider molecule formation, we will consider instead the 
second branch in Fig. ^ for the renormalization. The 
renormalization for positive scattering length will only 
be valid when two-body potential that does not support 
a bound state. The energies used for the renormaliza- 
tion are Eqm = 3/2fko and the energy branch where 
\/2hw < E re i < 5/2huj. This branch of solutions is a 
smooth curve which gives the correct non-interacting en- 
ergy at a — 0. 
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FIG. 1: (Color Online) Spectrum E re i as a function of ao- 
The dashed blue line corresponds to the energy branch se- 
lected for the renormalization procedure developed in this 
study. 



B. Mean- field solution 

Next we determine a renormalization function £(fc/ao) 
which, applied to the system of two particles in a trap 
and using the HF approximation, yields the exact results 
obtained in the previous section. 

To obtain the HF solution of two opposite-spin parti- 
cles in a trap, we utilize a product wavefunction having 
the same orbital for both particles. Thus, the two-body 
spatial wavefunction is 

*(ri,r 2 ) =V(ri)^(r 2 ) (6) 

and the spin part is antisymmetric. We introduce this 
trial wavefunction into the Hamiltonian and obtain the 
energy functional, 



Sty) = J (V(r) + l -mu 2 v 2 ^ ^(r) 

mkf ) 

Minimization of this energy functional determines the 
ground state energy and wavefunction. The minimiza- 
tion is done with res pec t to the orbital tp(r) as in a stan- 
dard HF procedure. 28J But prior to carrying out this 
minimization of Eq. J7J, we must choose how to eval- 
uate kf. Since its formal definition is kf = (6tt 2 p) 1 / 3 , 
this means that kf depends at each r- value on ^(r). For 
many-particle systems, we would use local density ap- 
proximation to evaluate kf. But the application of a 
local density approximation for a system of two particles 
does not seem physically correct, so, for two particles 
we consider the expectation value of kf, to be the more 
appropriate quantity. 

k f = J kf(r)^(r) 2 dr= /"(67r 2 ^(r) 2 ) 1/3 ^(r) 2 dr (8) 

The minimization procedure leads to a Schrodinger-type 
equation, where ip 2 ( T ) is the 1-particle density: 

(~—V 2 + -muj 2 r 2 + 
\ 2m 2 

STr^TT 2 ) 1 /^ / C(k f _ao)ao _ C(fc/aoA , (r)2/3 
3m [ kf k) ) 

+ ^f^(r) 2 W) = e*(r) (9) 
mkf J 

where e is a Lagrange multiplier which represents the 
chemical potential. The relation between e and the en- 
ergy is not as straightforward as in the HF case, owing 
primarily to the appearance of C(^f a o)- It should be un- 
derstood that ifj(r) 2 ^ 3 is supposed to be evaluated on a 
branch for which it is real and positive everywhere. Here 
and in the following, C'(x) = dC,{x)/dx. 

Equation JHJ corresponds to the GP equation for 2 
particles with a renormalized scattering length a e ff = 
C(kfa )/kf. After solving Eq.@ we use 10 to evaluate 
the energy. The basic idea is, for any chosen bare two- 
body scattering length ao, to find £(fc/ao) so that the 
energy of the ground state of 10 matches exactly the ap- 
propriate energy of fijl. From our numerical experience, 
the functional dependence of £ on fc/ao appears to be 
uniquely defined by the set of equations l@J [SJ and 
©■ 

There are two self-consistent procedures involved in 
this calculation. To solve Eq. © we follow the standard 
HF procedure, in which we adopt noninteracting solu- 
tions as the initial guess for the orbitals, after which we 
iterate Eq. © until convergence is achieved. For this 
procedure, we need the functional form of £(fc/ao) and 
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C'(fc/ao) over a range of kfao values since kf is chang- 
ing in each iteration. This means that we cannot find 
the exact renormalization function C(^/ a o) at any fixed 
value of ao without knowledge of the functional form of 
C(fcyao) at nearby values. To solve this problem we cal- 
culate C(fc/&o) self-consistently over the entire range in 
kfdo that is of interest. First we select a set of scattering 
length ao values which cover the entire range of interest. 
For an initial trial dependence (k / ao) we solve Eqs. 
© and J7J at each ao, obtaining the energy E, kf, the 
wavefunction and e. Then, to obtain a new (kfao), 
we look for the value of scattering length So for which 
Eexact^flo) = E, and we generate a new renormalization 
function that satisfies C '(kfao) = (^(kfdo). The mod- 
ification of the renormalization function is evidently in 
the abscissa rather than in the ordinate. This is a con- 
venient way to approach this calculation. Once we have 
carried out the matching procedure onto the whole set of 
scattering length ao values, we generate the next itera- 
tion for C (^»/ ffl o) and its derivative by interpolation. 

In the next iteration, ^°\kfao) is replaced by 
£W(fc/Oo) and we repeat the energy matching step for 
the whole set of ao-values. This procedure is repeated 
a few times until it converges to give a single correct 
renormalization function function £(fe/ao). Note that 
this iterative procedure determines a "numerically exact" 
renormalization function f(fc/ao). Because the iteration 
procedure is efficient, in 5 iterations we obtain 9 digits of 
agreement between E exact and E over the entire ao range. 
It is important to introduce a sensible initial trial renor- 
malization function (^{kjao). Many trial £(°)(fc/ao) 
functions, like Q°'(kfao) = kfao, would produce collapse 
of the two-fermion wavefunction for large and negative 
ao. To avoid this collapse, we propose an initial trial 
£(°)(fc/Oo) which is close to the correct £(fc/ao), this is 
done by choosing a qualitatively correct functional form 
with a few free parameters and we then find the set of 
parameters that best reproduce the exact two-body en- 
ergies. 

The final numerical results obtained for the renor- 
malization function C(fcyao) are accurately approximated 
by the monotonic functional form £o(kfao) = A + 
J5arctan(C + Dkfao), where A and B are chosen to 
have the corresponding maximum and minimum val- 
ues at ao — > ±oo and C and D are chosen to obey 
£(kf ao) — > kfao for kfao << 1- The maximum value 
is ( max = 2.182 and value is C, mm = -1.392, this leads 
to A = 0.395_and B = -1.138. To get the correct 
behavior for kfao << lj this m turn requires C = 
a,rctim(-A/B) w 0.362 and D = -(1 + C 2 )/B « -0.994. 
Thus there are only two independent parameters A, B 
to be specified at this level of approximation. Figure 
121 compares our numerical results for ^(k/ao) with this 
arctangent approximation, 

Co(k f a ) = 0.395-1.138 arctan(0.362-0.994fc / a ). (10) 

Figure |3 displays the fractional error in Co(kfao) defined 




kfag 

FIG. 2: Effective scattering length £(fc/ao) (circles) and its 
analytical approximation Co(kfdo) (full line) 



TABLE I: Exact numerical values of £(x). 
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«*) 


X 


C0=) 


x 




— oo 


-1.392 


-1.6582 


-0.82404 


4.9199 


1.9248 


-11.937 


-1.2818 


-1.3097 


-0.73716 


5.2872 


1.9426 


-11.394 


-1.2767 


-0.96693 


-0.62355 


5.6545 


1.9581 


-10.85 


-1.2712 


-0.63317 


-0.47254 


6.0219 


1.9717 


-10.306 


-1.2651 


-0.31351 


-0.27173 


6.3892 


1.9838 


-9.7626 


-1.2584 


-0.014329 


-0.014242 


6.7566 


1.9946 


-9.2193 


-1.251 


0.26087 


0.2863 


7.1239 


2.0042 


-8.6758 


-1.2428 


0.51684 


0.59122 


7.4913 


2.0129 


-8.1326 


-1.2336 


0.76339 


0.86101 


7.8586 


2.0208 


-7.5895 


-1.2231 


1.0072 


1.0787 


8.226 


2.028 


-7.0467 


-1.2113 


1.2507 


1.2474 


8.5933 


2.0346 


-6.504 


-1.1977 


1.6166 


1.4309 


8.9607 


2.0406 


-5.9617 


-1.182 


1.9829 


1.5585 


8.0056 


2.0238 


-5.4197 


-1.1635 


2.3497 


1.6508 


8.9117 


2.0398 


-4.8782 


-1.1416 


2.7166 


1.7201 


9.523 


2.04896 


-4.3373 


-1.1153 


3.0837 


1.7737 


9.8669 


2.0536 


-3.7973 


-1.0829 


3.4509 


1.8164 


10.895 


2.0657 


-3.2585 


-1.0423 


3.8181 


1.8512 


11.311 


2.0699 


-2.7216 


-0.99004 


4.1853 


1.88 


11.998 


2.0763 


-2.1875 


-0.92043 


4.5526 


1.9042 


+oo 


2.182 



as (C(fc/ao) — £o(fc/ a o))/C(^/ a o)) showing a maximum 
error of approximately 5%. 

Now that the renormalization function has been deter- 
mined, other observables can be tested for the two par- 
ticle system. Interestingly, there is a numerically exact 
agreement between the the external trap potential energy 
expectation values measured with the exact wavefunc- 
tion and with the mean-field renormalized wavefunction. 
However, the one-particle density profiles calculated us- 
ing the exact wavefunctions and the mean-field renor- 
malized wavefunction are only in qualitative agreement, 
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kf do 

FIG. 3: Fractional error in our analytical approximation to 
the numerical renormalization function, fo(fe/Oo)- 



e.g. for scattering lengths of large magnitude, where 
|fc/a | >> 1. 
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FIG. 4: Ratio of the total energy to the non-interacting 
energy, for a spherically-trapped two-component degenerate 
Fermi gas in the large N limit. The circles correspond to HF 
calculations for 2280 particles using £(fc/ao), while the solid 
line corresponds to either the variational solution, eq. or 
the TF solution, eq. (1251 , (the curves are indistinguishable on 
the scale of the figure) using the approximate renormalization 
function Co(kfdo)- 



III. APPLICATION TO MANY-PARTICLE 
SYSTEMS 

This section presents different many-particle approxi- 
mations for which the renormalized scattering length can 
be used. The renormalization procedure is designed to 
be used in the Hartree-Fock approximation, however we 
will see that simpler approximations like Thomas-Fermi 
(TF) or a variational trial wavefunction will yield equally 
effective results in the large N limit. The variational 
wavefunction we will use is the noninteracting wavefunc- 
tion rescaled in the radial direction by a factor A as in 
Ea. (|ll|l . where A is the variational parameter. As an ex- 
ample, we can see in Figure0]a comparison for the ground 
state energy of a two-component Fermi gas in an spheri- 
cal trap in the large N limit. The result obtained with the 
approximate Co(^/ao) in conjunction either with a vari- 
ational trial wavefunction or else with the TF method 
are in good agreement with the full HF calculation with 
the exact £(fc/ao). The difference between the results is 
mainly used by the replace of the exact £(fc/a ) by the 
approximate £o(fc/ a o)- If we use the exact £(fc/ao) for 
all the methods, the energies agree in at least 3 digits. 
In systems having a small number of particles, the HF 
method is, of course, most reliable. 



A. Variational 

The simplest approximation |29j utilizes a trial wave- 
function that is a simple radial rescaling of the noninter- 



acting wavefunction: 

*A(ri,r 2 , ...,rjv) = 3jy/2 ^ jvj(r x /A, r 2 /A, r N / A) 

(11) 

The expectation value of the renormalized Hamiltonian 
can be separated into two terms, E(X) = EhqW + 
E in t(X,a ), where 

i ^ ' 

E int (X,a ) = (* A | V 47Th2aeff S(r z - ry) |* A ) . (12) 
* — ' m 

i<_i' 

The energy of this trial wavefunction is calculated 
as a function of the variational scale parameter A for 
the renormalized Hamiltonian The non- interacting 
wavefunction is a Slater determinant formed with the oc- 
cupied spin-orbitals. The Eho is simple to calculate, as 
it requires only a change of variables to determine the 
A-dependence in Eq. 1121 in conjunction with the known 
results of the non-interacting ground state. 

E HO (\) = E NI (J^ + (13) 

The interaction energy Ei nt can be written in the follow- 
ing form, when we apply the renormalization locally as a 
function of the density: 

E int (X,a ) = —J fc A (r) plWr (14) 
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In this equation p is the density of one spin component 
and kf(r) = (67r 2 p^(r)) 1 / 3 . In the large N limit, the den- 
sity of the non-interacting wavefunction can be replaced 
by the TF density of the noninteracting system [30( . The 
density corresponding to our trial wavefunction is a sim- 
ple radial rescaling, whereby the density in the high-N 
limit is: 



/67V 



px(r) = { ^ 2 <o^ 




2a£ o A 2 (3JV) 1 /3 



3/2 



if r 2 < R 2 

otherwise 
(15) 

Here N is the total number of particles and R c = 
V / 2ah A(3A r ) 1 / 6 is the radius of the Fermi gas. In the 
large N limit, the total energy can be expressed in units 
of the noninteracting energy, 



A 2 
2 



k° f a S 



(16) 



Here A is the scaling parameter, = v / 2(3A r ) 1 / 6 /a/j is 
the Fermi momentum of the non-interacting system at 
the trap center, and F is 

F(7) = ^ / " x2)5 ' 2x ^ ( 7v/l = x2 ) dx - (17) 

This function must be calculated numerically unless fur- 
ther approximations are made. The energy results ob- 
tained using Eq. (|lu|) are shown in Figure In the 
unitarity limit, the behavior can be calculated exactly: 



F( 7 



-oo 



4 4 C r ' 



9tt 2 



(1 - x 2 f' 2 x 2 dx 



5C 



9tt 



(18) 

Three curves predicted by Eq. (fTT?f> are shown in Fig. 
For the entire range of interactions, the energy of the sys- 
tem (the minimum of the curve) remains finite, ranging 
from 0.713£jvj to 1.33-E/vj, This shows how our renor- 
malization circumvents the collapse that would occur for 
the bare Fermi pscudopotcntial. 



B. Thomas- Fermi results 

In this subsection we will review the TF approxima- 
tion, using the renormalization function. The TF ap- 
proximation has been used to study a 2-component Fermi 
gas with zero-range pseudo-potentials |3l] , but no renor- 
malization has been considered. 

Thomas-Fermi is of course a local density approxima- 
tion. At each position inside the trap, the wavefunction 
is approximated by a Slater determinant of a set of plane 
wave orbitals, i.e., the orbitals are characterized by four 
quantum numbers, the vector momentum k and the spin. 
The orbitals are filled uniformly up to a level fc/(r) which 
is the same for spin up and spin down fermions. The 
value of kf(r) will depend on the distance r from the 




FIG. 5: Total energy as a function of A, in units of the 
noninteracting total energy. The solid curve corresponds to 
fetdo = 0, the dashed curve corresponds to k^ao = — oo, 
and the dotted-dashed curve corresponds to k ( jao — oo. The 
minimum of the energy functional for k'jao ~ — oo occurs at 
A = 0.844, which represents the ratio between the cloud ra- 
dius at unitarity and the noninteracting cloud radius. 



trap center, and on the number of particles in the sys- 
tem. For a uniform system, the value of kf is a constant 
that characterizes the density of the system. 

To calculate the local energy we need to sum over all 
the states at that position. For example, the kinetic en- 
ergy term K for one-spin component is: 



vh 2 



2to(2tt) 



20tt 2 



rn 



(19) 



Here V is the volume of integration which will disap- 
pear when we consider the local energy. This volume is 
small in comparison with the external potential (in this 
case the trap) characteristic length but is big enough to 
contain many particle. So, kf and V ex t can be consid- 
ered constant during the integration. The calculation of 
the expectation value of an external trapping potential 
is then straight forward and we obtain 



(V ex t) — -^r^Vextkj:, 



(20) 



which is just the number of particles times the exter- 
nal potential at that position. In the case of the inter- 
action term, the two different spin components can be 
calculated by considering the particles as indistinguish- 
able, and calculating the direct and exchange terms. But 
since the interaction is a delta function, we can consider 
the particles as distinguishable with interactions only be- 
tween different species and obtain the same result. This 
latter procedure is easier and we only need to consider 
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the direct term. 

(Vint) =^(kk'|^„ t |kk') 



kk' 

1 eff 



1 Ana pff h 2 rkf rk 



(27r) e 



in 



f r K f 
o Ja 



5(-K-x')dx'dxd 3 k'd 3 k 



v ^a eff h 2 kl ( fc/ )3 = y 4ntf fc /C( fc /»o) (21) 



67T 2 67T 2 



TO (67T 2 ) 



2\2 



For the case of two equally-numerous spin components 
the local energy (per unit volume) is 



,. kl 4-irh 2 k 5 f((k f a 

e (kt) = E {kl)/ v = -^ + v^ : 



it k) 



TO (67T 2 ) 2 

In an infinite uniform system, where V ex t — 0, the energy 
is: 



^_fc|_ 47rfi 2 k 5 f C(k f a ) 
2m 57r 2 to 367T 4 



(23) 



The ratio between the total energy and the non- 
interacting energy has a simple form and only depends 
on fc/do, 

S(k f )/6Mk f ) = l + 1 -^^. (24) 

Using Eq. 1221 we can construct our energy functional by 
integrating the local energy over all space. 



E 



dr ( £±+- + V ext {r)-- 



2m 57T 2 



3vr 2 

Anh 2 k f {r) 5 C(k f a ) 

TO 367T 4 



(25) 



To find the ground state we have to minimize the en- 
ergy under the constraint that the number of particles is 
fixed. This constraint can be implemented by introduc- 
ing a Lagrange multiplier /io, usually called the chemical 
potential. So, the minimization of Eq. I|25|l for fixed 
number of particles is reduced to the minimization of 

A = E- l i N = E-fi J dr^^, (26) 

where variational parameter is kt (r) . The necessary but 
not sufficient condition for fe/(r) to minimize A is that 



dk 
dkj(v) 



= 0. 



(27) 



This condition leads to a relationship between the local 
chemical potential, defined as £t(r) = (1q — V ex t(v), and 
the local Fermi momentum fc/(r), 



h 2 k 2 (r) f 10 

air) = — -r 1 + 7— 

PW 2m V 9tt 



,(fc / (r) ao ) + ^|^C , (fc/(r)ao) 

(28) 



The value of po fixes the number of particles and, with 
this relationship, we can calculate the density profile and 
the energy of the system. Figure El shows the chemical 
potential dependence on fc/ao obtained with the renor- 
malization function and with other models. In Figure 0] 
we can see the energy obtained using eqs. <|25I28|) in the 
large N limit. 

In the unitarity limit, oq — > — oo, we obtain 



h 2 k 2 ( 10 
2m \ 9ir 



(29) 



At unitarity, when the scattering length is much larger 
than the inter-particle distance, the only relevant param- 
eter is the density |3^,|3^|. Dimensional analysis suggests 
that \i oc p 2 ' 3 cx k 2 . The expected relation between fi and 
kf is usually written as 



21.2 



K l k 



2to 



£ (l+/5). 



(30) 



From our calculations this relations appears naturally 
with a coefficient (3, which is an universal parameter, of 
{3 = 10C mm /97r = -0.492. This parameter f3 has been 
studied from many different perspectives. Table ITT1 shows 
different experimental and theoretical values of f3. Inter- 
estingly, our value is consistent with most experiments. 
To measure (3 experimentally, Hulet and Thomas groups 
measure the size of the cloud and compares it with the 
noninteracting cloud. For, example, the result obtained 
by Hulet group is Rjj/Rni — 0.825 ± 0.02 which com- 
pares well with ours Ru/Rni — 0.844 (Fig. 0, obtained 
using both variational or TF calculations we obtain. 




-2 -1 

l/kfCLQ 



FIG. 6: (Color Online) Chemical potential in units of the 
Fermi energy. The black solid line is the prediction obtained 
with the renormalization function. The black circles repre- 
sent the BCS prediction and the blue dashed curve is the 
prediction obtained in Ref. .34] 



It is well established |33j, l42t |43j that an ultracold two- 
component Fermi system exhibits superfluidity. Even 
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TABLE II: Experimental and theoretical predictions of f). 



p 


Experiments 


-0.54(5)" 


-0.64(15)'' 


-0.49(4) d 


Experiments 


-0.68±°;^ 


-0.54± ; il e 




QMC 


-0.58(1/ 


-0.56(1) 9 




Pade asymptotes 


-0.674 ft 


-0.432 i 




Green's function 


-0.545 J 


-0.599 fc 




BCS 


-0.41 ! 






Other methods 


-0.3 m 


-0.564" 


-0.492° 



a Ref. 
6 Ref. 
c Ref. 
d Ref. 
e Ref. 
•fRef. 
°Ref. 
h Refs 
'Ref. 
3 Ref. 
k Ref. 

'This is a well known result, sec for example Ref. 39 1 
m Ref.|4i| 
"Ref.0 

"Result obtained in this paper with the renormalization function. 



though our renormalization scheme does not explicitly 
consider superfluidity, it reproduces a number of prop- 
erties of the Fermi gas sensibly, including the equation 
of state and the chemical potential. Consequently, these 
results can be used in the hydrodynamic theory to ex- 
tract information about dynamics of the system, like the 
speed of sound or normal modes of excitation. For ex- 
ample, the speed of sound in a uniform two-component 
system is given by (3^ 



h d 
m dp 



,dE/N\ 
dp J 



(31) 



Using Eq. (|24[) we can thus evaluate the speed of sound, 
which generates the results shown in Figure[7] The speed 
of sound results reproduce the expected limiting behav- 
iors. In the noninteracting limit v — Vf/y/3, while at 
unitarity v — Vf^J (1 + /3)/3 ^ij. This is one example 
of a nontrivial observable quantity for this system that 
can be predicted by this renormalization technique. A 
comprehensive study of other observables based on this 
approach will be left for future publications. 



C. Hartree-Fock method 

The HF method for a many-particle system is an ex- 
tension of the two particle calculation done in Section 
II. The variational parameters are the orbitals and the 



0.55 



0.5 



0.45 



0.4 



-10 



-8 
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FIG. 7: The speed of sound is shown in units of the Fermi 
velocity Vf = hkf/m for an uniform two-component Fermi 
ffclS, £LS a function of fc/ao- 



energy functional is: 



( N/2 ( fi 2 1 \ 

2$>(r) (^-|-V 2 + i™c 2 r 2 J^(r) 



Sty) 



p(rf)dr. (32) 



4TTh 2 k f a eff (k f (r)a ) 
kf{v)m 



In this approximation, the renormalization is done lo- 
cally, kf(r) — (67r 2 p(r)) 1 / 3 . The minimization procedure 
is lengthly and straight forward, so it will not be pre- 
sented here. By minimizing with respect to the set of 
orbitals ipi, we obtain a set of nonlinear HF equations. 
These are solved self-consistently. Figure 0] shows results 
for the HF energy of 2280 particles. This approxima- 
tion is particulary useful for systems with small number 
of particles, for which the TF approximation has limited 
applicability. In Sec. IV below, this method is used to 
obtain the energies of 8 fcrmions in a trap. 



IV. RESULTS 

To compare the predictions based on our renormalized 
scattering length with other methods we have carried out 
fixed node diffusion Monte Carlo (FNDMC) simulations 
for equal mixture of different-spin fermions. Interac- 
tions are considered only between different-spin fermions, 
which are treated here as distinguishable particles. The 
interaction potential is a purely attractive gaussian and 
its width d is chosen so that pd 3 « 10~ 4 . The method 
follows closely the approach of j3|| . 

The ground state wavefunction of two particles in a 
trap can be separated exactly into a Jastrow term and 
non-interacting orbitals 



* P (ri,r 2 ) = ip(ri)^(r 2 )J{r 1 - r 2 ), 



(33) 
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where ip is the non-interacting ground state orbitals. 
We will show that this relation is valid when the center 
of mass wavefunction and the noninteracting orbitals 
are gaussians, like in the trapped. For this discussion 
H = m = u) = 1 and we will not consider normalization 
factors. If wc define Rcm = (l/2)(ri + r 2 ) and 
r = (ri — r 2 ) we know that the pair wave function 
can be separated in center of mass and relative co- 
ordinate term, \I/p(ri,r 2 ) = ^(Rcm)^! - ). Also, we 
know that the center of mass wavefunction is unaf- 
fected by the two-body interaction and is $(Rcm) = 
exp(-R 2 CM ) = cxp('-r 2 /2)cxp(-r 2 /2)exp(r 2 /4) = 
-0(ri)"0( r 2) exp(r 2 /4), so the total wavefunction 
is *(R C Af)0(r) = V(i'i)V'(i'2)exp(r 2 /4)(/ 1 (r) = 
-0(ri)-0(r 2 )J(r), where J(r) = exp(r 2 /4)0(r). The 
evaluation of the relative coordinate wavefunction, (j)(r), 
requires in general of a numerical calculation. 

We use this Jastrow term to construct the many- 
body wavefunction. This wavefunction is usually called 
Jastrow-Slater wavefunction. Here, i and i' correspond 
to different-spin fermions. The non-interacting wavefunc- 
tion is a Slater determinant formed with the harmonic 
oscillator orbitals. 

*(ri,r 2 , ...,ry) = J(?i r i ,)-^ NI (r 1 ,r 2 , ...,r N ) 

a' 

We have also used a BCS-type many-body wavefunction 
constructed with pair wavefunctions for FNDMC 

*Oi,ri, -,r N/2 >) = 

A {* P (ri, ri-)* p (r 2 , r 2 / )...*„( r JV/2 7 r N/2' )}, (35) 

where A is the antisymmetrizer operator. This trial 
wavefunction leads to good results on the BEC limit but 
in the BCS regime the Jastrow-Slater wavefunction pro- 
duces lower energies. 

We have calculated FNDMC energies for 8 particles in 
the BCS side of the crossover. In Fig. [SJthese energies are 
compared with HF calculations including the first- and 
second-order corrections in the k/ao expansion (2?|], and 
then also with full HF calculations using the renormalized 
scattering length directly. 



E mt /N 



h 2 kj fkfdo 6(11- 2 In 2) 



3tt 



1057T 2 



(Mo) 



(36) 

The idea of using this type of expansion to construct 
energy functionals has been applied for bosons 0, . 
The expansion l|36(l can be introduced locally in varia- 
tional treatments, which yields an energy functional, 



£{tl>) 



/ N/2 
] ( 2 5>(1 



H 2 I 

V 2 + -muj 2 v 2 I ibJr) 

2m 2 ' 



+ 4^ + a , 12(ll-21og( 2)) /3p(r)7/3 ) dF 

m 1057r z 




FIG. 8: (Color Online) The total energy of 8 fermions in a trap 
is shown in oscillator units as a function of k° f a . FNDMC 
results are shown in open red circles while full blue line cor- 
responds to HF results using the effective renormalized scat- 
tering length. The dashed and the dash-dotted curves cor- 
respond respectively to solution use first-order, or first- and 
second-order terms, in an expansion into powers of fc/ao- See 
the discussion of Eqs. 1361371 . 



where p = p| = pj. If we only consider the first term in 
Eq. 1361 we obtain the Fermi pseudo-potential contribu- 
tion. 

To study the weak interacting limit, previous authors 
pTol l31| have considered the Fermi pseudopotential ap- 
proximation, which is only the first term in the energy 
expansion ll.'Uill . Applying the expansion (|36|l in the local 
density approximation is a convenient way to introduce 
higher order corrections to mean field theories. We can 
obtain an expansion of the density-dependent renormal- 
ization function using Eq. I|36|) . in this case 



C(fc/ao) = kfa + 



6(11 - 2 In 2) 
35^ 



{k s a Q f 



(38) 



(37) 



Insertion of this result into Eq. (JIJ , with the local density 
approximation and a Slater determinant wavefunction, 
gives Eq. i|371l . 

A power expansion of £(fc/ao) obtained by the renor- 
malization method should agree with this expansion. The 
first term is reproduced exactly but the second one is 
only in qualitative agreement. While the coefficient of 
the second order expansion in Eq. (|38|l is approximately 
0.525, in the density renormalization from Section II the 
coefficient is 0.422. This disagreement may be due to 
the level of approximation of the density renormalization 
procedure. 

We find very good agreement between the mean-field 
results calculated using the renormalized interaction de- 
veloped in this paper, and the FNDMC (Fig. |SJ). The 
variational methods including the perturbative correc- 
tions (Eq. I37fl show good agreement in the small fc/ao 
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region, deviating from the FNDNC results when the cor- 
rections to the expansion (|36[) become important. 



1 




-6 -5 -4 -3 -2 -1 

1/fc/Oo 

FIG. 9: E/Eni for an infinite homogeneous Fermi gas fc/ao 
in the mean-field approach (solid line). The dashed curve 
corresponds to the local density BCS solution and the circles 
correspond to FNDMC results obtained in [3(J. 

It is also possible to make comparisons with other 
quantum Monte Carlo calculations. Astrakharchik and 
coworkers |36| have studied a homogenous doubly degen- 
erate Fermi gas using FNDMC methods. In their calcu- 
lations, they considered up to 60 particles. We compare 
the energy of this system as obtained using the density 
renormalization procedure in Eq. 1)24(1. A comparison be- 
tween the two calculations and local density BCS result 
USE! is shown in Fig. J3J. 

V. CONCLUSIONS 

It is the goal of most many-body theoretical studies 
to derive predictive power for numerous observables of 
interest, using simpler methods that bypass the actual 
calculation of this "true" ground state wavefunction for 
the trapped atomic gas. At the heart of many such treat- 
ments are the following two steps: (i) replacement of the 
two-body potential energy by a zero-range Fermi pseu- 
dopotential, followed by (ii) a mean-field wavefunction 
ansatz and the computation of observables. The basic 
level of description for a Bose gas incorporates no cor- 
relations whatsoever. For a Fermi gas, correlations are 
generally treated at either the bare minimalist level of ex- 



change correlations alone, using a single Slater determi- 
nantal wavefunction. A more sophisticated level is often 
considered for a system of mutually attractive fermions, 
which are frequently described with BCS-type correla- 
tions built into the description. One way of visualiz- 
ing the value of a Fermi-type zero-range pseudopotential 
adopted in most such theories is to remember that it 
has been specifically designed to give a meaningful in- 
teraction energy for each pair of particles even when the 
wavefunction structure is too simplistic to incorporate 
any appreciable correlations. 

The present article presented an alternative implemen- 
tation of this general philosophy. We developed a pro- 
cedure for renormalizing the coefficient of a zero-range 
potential, based entirely on an analysis of the nonper- 
turbative two-body system solved first with and then 
without wavefunction correlations. When we applied this 
procedure to the many-body Fermi gas, it gives agree- 
ment with the standard dilute gas limit, an important 
prerequisite for any realistic theory. But in addition, it 
is able to treat higher densities n, including the regime 
|nao 3 | >> 1. We studied a number of observables that 
have been explored both experimentally and theoretically 
in the BCS-BEC crossover regime, and found good agree- 
ment using our renormalized Hartree-Fock approach all 
the way to the unitarity limit ao — > oo. Perhaps surpris- 
ingly, this good agreement is achieved without incorpo- 
rating explicit BCS-type correlations into the many-body 
wavefunction. One result of this study is an approximate 
expression for the renormalization function in closed an- 
alytical form that may prove to be useful in other studies 
of the two-component degenerate Fermi gas. Another in- 
teresting result is that at unitarity, the chemical potential 
exhibits the expected density dependence characterized 
by the parameter (3 = —0.492, which, interestingly, is 
consistent with recent experiments 0, S S 3 ■ 

In order to study the complete BCS-BEC crossover, 
future improvements of this theory should include a more 
flexible many-body wavefunction which can represent a 
Fermi gas in the weak interacting region and a gas of 
Bose molecules in the BEC region. 
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